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i.  INTRODUCTION 


1 . 1  PURPOSE,  SCOPE,  AND  LIMITATIONS 

This  is  Volume  X,  Lunar  Gravity  Analysis,  of  a  documentation  series  for 
the  TRACE66  Trajectory  Analyses  and  Orbit  Determination  Program.  In 
particular,  it  is  intended  as  a  technical  reference  for  the  Lunar  Gravity  Field 
Analyzer  of  TRACE66.  Significant  analysis  procedures  and  techniques  used 
to  determine  lunar  gravitational  constants  and  reconstruct  lunar  satellite 
orbits  are  described. 

The  symbols  and  notations  of  coordinate  system  transformations  associated 
with  reference  frames  and  equations  throughout  the  document  are  defined 
and  discussed  in  Section  2.  Equations  of  motion  and  variational  equations  are 
expressed  in  Section  3,  with  emphasis  on  the  description  of  each  model  used 
to  evaluate  the  various  accelerations  acting  on  the  vehicle. 

In  Section  4,  the  observation  measurement  model  for  long-count  two-way 
doppler  data  is  described.  Sections  5  and  6  present,  respectively,  the 
details  associated  with  the  accumulation  of  the  normal  matrix,  and  the 
solution  of  the  normal  equations. 

Techniques  used  to  obtain  a  composite  lunar  gravity  field  are  described 
in  Section  7,  while  features  peculiar  to  the  TRACE66  Lunar  Gravity  Field 
Analyzer  are  discussed  in  Section  8. 


-1- 


1.2 


HISTORICAL  BACKGROUND 


The  fundamental  purpose  of  TRACE  is  to  simulate  orbital  motion  and  tracking 
operations.  In  the  past,  lunar  orbit  determination  was  performed  at  The 
Aerospace  Corporation  by  the  TRACE- DL  Orbit  Determination  Program 
(Ref.  1).  The  current  TRACE66  program  has  significant  capabilities 
that  surpass  those  of  TRACE-DL. 

I-3  TECHNICAL  BACKGROUND 

In  particular,  this  document  summarizes  the  software  developed  for  the 
dynamical  dete  -mination  of  lunar  gravitational  constants  and  reconstruction 
of  lunar  satellite  orbits.  Recent  analyses  and  investigations  (Refs.  2  and  3) 
of  NAoA  Deep  Space  Network  doppler  observations  of  Lunar  Orbiters  have 
utilized  TRACE66  and  auxiliary  software  for  dynamical  fitting  and  composite 
model  determination. 

1. 4  SUMMARY 

This  document  describes  both  the  TRACE66  Lunar  Gravity  Field  Analyzer 
and  external  auxiliary  software  that  determine  lunar  gravitr tional  constants 
and  reconstruct  lunar  satellite  orbits  from  NASA  Deep  Space  Network  2~way 
doppler  data.  The  force  model  used  to  generate  lunar  trajectories  is  a 
combination  of  perturbative  accelerations  resulting  from 

•  Gravitational  potential  due  to  a  surface  distribution 
grid  (7  50  points)  of  disk  or  point  masses  {p^,  0^,  M,  a^} 

•  Lunar  gravitational  potential  {p,  ^nmi  ^nrr-} 

•  Gravitational  potential  due  to  planetary  bodies 
(e.g.,  Sun  and  Earth). 

Principal  features  of  interest  here  are  the  following  capabilities: 

•  Mass  parameter  selection  by  orbital  arc 

•  Doppler  residual  generation  with  printer  plot  versus 
vehicle  seienographic  ground  track 
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•  Residual  and  edit  summaries 

•  Self-contained  orbit  determination  function  if  the 
number  of  parameters  does  not  exceed  100 

•  Normal  matrix  (256  X  256)  accumulation  per  arc 
(exterior  to  TRACE66,  these  normal  matrices 
are  merg-.’d  and  the  resulting  linear  system,  of 
the  order  2000,  is  solved). 
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2.  SYMBOLS  AND  NOTATION 


The  common  symbols  used  throughout  this  document  are  presented  next, 
followed  by  a  discussion  of  the  basic  notation  associated  with  coordinate 
system  transformations. 


2.1  SYMBOLS 

•  time  differentiation  and  vector  dot  product  operation 

••  time  differentiation 

as  an  underscore,  indicates  a  vector;  for  example,  _r  is 
the  vehicle  position  vector 

A  increment  or  difference 


X  vector  cross-product  operation  and  product  sign 

absolute  value,  or  vector  magnitude  operation;  for  example, 
r  =  |r  j 

V  vector  gradient  operator 

{  )  indicates  a  functional  relationship  or  dependence,  such  as 

the  vehicle  acceleration  vector  £(£*£.*  y),  which  is  a  function 
or  £,£,  and  y.  Also,  used  to  denote  row  vectors;  for  example, 
r  =  (x,  /,  z) 

{  }  indicates  a  set  of  elements;  for  example,  {V,p,$,  X-} 

[  ]  a  matrix;  for  example,  r  =  [x,y,z]  (row  vector)  or 


r 


(column  vector) 


tntiM 
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n  o  01 


Superscripts 


i. ) 
c 
T 
n 


I  = 


0  i  0 
0  0  i 


3X3  identity  matrix 


particular  index  relative  to  a  set  or  vector 
computed  value 
matrix  transpose  symbol 
iteration  index 


Subscripts 
0,  1,2,3 
i»  j,k,« 

M 

c,  me 
T,  C 
<h 
nm 
m,  n 
x,y,  z 


index  denoting  forces,  components,  or  axis  of  rotation 

particular  index  denoting  members  of  a  set  or  components 
of  a  vector.  In  general,  i,  j , k,  £  =1,2,  .  .  . 

related  to  the  Moon 

computed  value,  and  measured  minus  computed  value 

transmitter  and  count  time 

local  horizon  reference  frame 

degree  n  and  order  rn 

denotes  index 

coordinate  axes 
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2.2 


COORDINATE  SYSTEM  TRANSFORMATIONS 


Lunar  gravity  analysis  requires  the  consideration  of  several  reference 
coordinate  systems  (see  Volume  II).  The  basic  transformations  of  interest 
here  relate  one  system  to  another;  they  are  described  by  the  following 
examples: 

In  the  Moon-Centered  Inertial  (MCI)  relationship  to  the  selenographic  or 
Moon  Fixed  (MF)  frame 


MF 

-TEE 


=  [M] 


MCI 

^■MEE 


where 

IvlOX 

— MEE  ~  position  vector  of  the  vehicle  referenced 

''"W;  J  to  a  mean  equator  and  equinox  of  some  base  date 

[M]  =  time -dependent  transformation  matrix  from  a 

mean  rectangular  inertial  frame  of  base  date  to  a  true 
selenographic  (moon  fixed)  frame 

[M]  =[S][P] 

=  [SP] 

[P]  =  time-dependent  precession  matrix,  which  transforms 

from  mean  equator  and  equinox  of  some  base  date 
to  mean  equator  and  equinox  of  date  (see  Figures 
1  and  2). 

[S]  =  time-dependent  transformation  matrix  from  mean 

rectangular  inertial  of  date  to  true  selenographic 
of  date 

—TEE  =  position  vector  of  the  vehicle  referenced  to  a 
true  equator  and  equinox  of  date. 
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In  the  selenographic  or  Moon-Fixed  (MF)  relationship  to  the  Moon-Centered 
Inertial  (MCI) 


-MCI 

-MEE 


=  [M'1] 


-MF 

-TEE 


where 

[M]  =  [M]  =  0  (by  assumption) 


..MF  ,  ,  , 

—TEE  =  acceleration  vector  of  the  vehicle  referenced  to  a 
true  equator  and  equinox  of  date 

[M]  *  =  time -dependent  transformation  matrix  from  a  true 
selenographic  frame  to  a  mean  rectangular  inertial 
frame;  note  that 


[M-1]  =  [SP]"1 


..MCI 

-MEE 


=  MCI  acceleration  vector  of  the  vehicle  referenced  to  a 
mean  equator  and  equinox  of  some  base  date  (this  is  the 
coordinate  and  time-keeping  system  in  the  chosen 
integration  system). 


In  the  Moon-Centered  Inertial  (MCI)  relationship  to  the 
Inertial  (ECI): 


Earth- Centered 


Figure  1.  The  Mean  and  True  Equator  and  Equinox  of  Date 


z 


OF  DATE 


Figure  l.  The  Moon- Fixed  System  and  the  True 
Equator  and  Equinox  of  Date 
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3.  EQUATIONS  OF  MOTION  AND  VARIATIONAL  EQUATIONS 


In  the  TRACE66  Lunar  Analyzer,  the  vector  equation  of  motion  is 

2 

»  V  ^ 

1  -  2-4  Li  (1) 

i=0 

where 

£  =  total  acceleration  of  vehicle 

£q  =  acceleration  due  to  discrete  masses  imbedded  in  the  Moon 

£  j  =  gravitational  acceleration  due  to  the  Moon 

£2  =  gravitational  acceleration  due  to  the  Sun  and  Earth 

Equation  (1)  is  numerically  integrated  in  a  Cowell  formulation  with  an  8th- 
order  predictor/corrector  Gauss- Jackson  differencing  method. 

Each  acceleration  considered  is  evaluated  in  its  appropriate  reference  frame. 
If  necessary,  it  is  rotated  and  accumulated  in  a  moon-centered  inertial 
(MCI)  frame,  which  is  then  used  as  the  integration  coordinate  system  (see 
Volume  II). 

The  variational  equation  for  any  equation-of-motion  parameter  y  (Ref.  3) 
is  expressed  by  the  vector  differential  equation 


I 


l 

I 


M 

9r 

or 

9r 

9r 

. 

^  ± 

9y 

Dr 

L  - 

dy 

,9r  . 

8r  3£(r,£,  y) 
dy  +  9y 


(2) 


Partial  derivative  matrices  9r/9£and  9r/9£  are  obtained  by  differentiating 
the  total  acceleration  with  respect  to  £  and  £.  The  first  two  terms  of  the 
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three  in  Eq.  (2)  account  for  the  implicit  appearance  of  y  through  the  functional 
dependence  of  jr  and  _r .  The  last  term  is  referred  to  as  the  nonhomo geneous 
term,  and  accounts  for  the  explicit  appearance  of  y  in  the  total  acceleration. 
Because  y  is  independent  of  time  and  the  derivatives  are  continuous,  the  order 
of  differentiation  in  Eq.  (2)  may  be  interchanged.  A  double  integration  with 
respect  to  time  then  yields  the  trajectory  partial  derivatives  8r_/3y. 

In  the  following  three  sections,  the  models  used  in  evaluating  the  accelerations 
of  Eq.  (1)  are  presented  in  detail;  also,  the  partial  derivative  matrices 
[9r/ 9rJ  and  [ 9r /  9rJ,  and  nonhomogeneous  terms  for  the  variational  equations, 
are  presented  where  applicable.  In  Section  3.4,  the  use  of  partial  derivatives 
is  discussed. 
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3.  1  ACCELERATION  DUE  TO  DISCRETE  MASSES 

In  TRACE66,  a  discrete  mass  force  model  can  be  represented  by  a  point  or 
disk  mass  formulation. 

3.1.1  Point  Mass  Acceleration 

Acceleration  due  to  point  masses  imbedded  in  the  Moon  is  given  by 


*0 =  771 


|£- 


(3) 


where 


th 


P.  =  the  selenographic  cartesian  position  vector  of  the  i 
point  mass 

th 

p^  =  the  mass  ratio  of  the  i  point  relative  to  the  Moon 
p  =  the  gravitational  constant  of  the  Moon 

[m]  s  the  transformation  matrix  from  the  mean  selenocentric 
frame  of  base  date  to  the  true  selenographic  frame 

_r  =  the  selenocentric  (MCI)  position  vector  of  the  vehicle 
where  £  =  (x,  y ,  z) 

th 

The  selenographic  position  vector  P.  of  the  i  point  mass  is  determined  by 

the  triplet  (p . .  X.),  where  p.  is  the  selenocentric  distance  of  the  point  and 

r  ‘  i  i  i  l 

<^,  X^  are  the  selenographic  latitude  and  longitude  of  the  point  mass.  In 

particular, 


P.  =  (p.  cos  $.  cos  X.,  p.  cos  $.  sin  X.,  p.  sin  $.) 

i  i  i  1  i  i  \  i  i 


(4) 
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The  mass  ratio  pu  is  the  differentially  correctable  point  mass  parameter.  The 
variational  equations  for  p..,  p^,  X^,  have  the  nonhomogeneous  terms 

siL  x- 

9|1i  h  |r  -  [M_1J  P.|3 


pp.[M-1]  8P.  f  .  8P.1  1-CM"1]^: 

'  3r  ■ 1 M' ,£i) '  ^rj  i£-Tm~]Tf 

p.p  .[  M_1]  3P  3P  1  r-[M_1]P. 

r-tL"1]?.!3  ^'3L(-'LM  ]^)“S^J]r.[M-1]Pl|Z 


pp.tM"1]  ap. 

r  „-lln  !3 

-  [M  JR j  1 


71— 3  ~Sk7  "  3  (-  -  [^"A]jPi)  *  -axT-  - - j— 

]  P.  i  L  1 J  |£  M 


'P.]r-  [M"1]  R 
\\  |r  -l  M"1]  P.|2 


where 


ap.  T 

=  (cos  $.  cos  X.,  cos  $.  sin  X..,  sin  $. ) 
00.  1  1  11  i 


ap.  T 

=  (-pi  sin  $.  cos  X^,  -pi  sin  ^  sin  X^  pi  cos  &) 


9P.  T 

-  (-p.  cos  $.  sin  X.,  p.  cos  $.  cos  X.,  0) 

0\,  rl  1  11  1  1 


-14- 


The  derivatives  of  the  point  mass  acceleration  with  respect  to  the  vehicle 
position  _r  and  velocity  _r  are  given  by  the  matrices 


=  -H- 


y _ h _ 


-  [M"1]Pi)T 
|r  -  [M'J]P.|2 


(6) 


and 


=  [Oj 


(7) 


where  1  =  3X3  identity  matrix. 

3.1.2  Disk  Mass  Acceleration 

Acceleration  due  to  disk  masses  is  formulated  as  one  resulting  from  a 
limiting  form  cf  zero-thickness  disks  (Refs.  4  and  5). 

The  geometry  used  to  relate  the  disk  masses  to  the  lunar  center  is  given  in 
Figure  3. 

The  disks  are  considered  to  be  infinitesimally  thin,  planar,  and  usually 
located  near  the  surface  of  the  Moon,  the  mean  equatorial  radius  of  which 
is  a^.  Now  by  letting  n  be  normal  to  the  sphere,  coinciding  with  the  axis 
of  the  disk,  the  local  horizon  system  z has  its  origin  at  the 

center  of  the  disk  where:  z ^  is  parallel  to  n;  x^  is  directed  south  on  the 
plane  normal  to  n;  and  *s  east  on  the  same  plane.  The  selenographic 
system  (x,  y,  z)  origin  is  located  at  the  center  of  the  Moon,  where  z  is  along 
the  polar  axis,  and  x  is  in  the  direction  of  the  intersection  of  the  equatorial 
plane  and  the  zero  meridian. 
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The  position  and  acceleration  vectors  in  the  local  horizon  system  are 
expressed  as 


-lh  =  (xfh’  yfh’  zih)T 
-ih  =  (xfh'  yjf h*  Z^h)T 


(8) 


The  selenographic  position  vector  f\,  of  the  i^  disk  mass  is  the  same  as 

Eq.  (4).  Therefore,  the  transformation  of  an  arbitrary  selenographic  vector 

til 

— MF  *°  ^l6  *oca^  horizon  frame  specified  by  the  i  disk  mass  is  given  by 

-l  h  =  R2  (z  ’  ^)R3(K)-MF  "  ^  (9^ 


where 


Ml  -  *)  = 


R3!M  - 


sin  0 

0  -cos 

0 

0 

.1  0 

COS  0 

m 

0  sin 

0 

» 

cos  X 

sin  X 

0 

-sin  X 

cos  X 

0 

0 

0 

1 

0 


P 


] 


T 


and  0,  X  are  the  latitude  and  longitude  of  the  mass  disk,  respectively. 
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/ 


Figure  3.  Disk 


-1 


The  components  of  the  acceleration  vector  due  to  the  mass  disk  in  the  local 
horizon  system  arc 


£h 

th 


*  * 

y lh  3  ji 
y2h  "  2  a3 


ra/iT 

a2+K 


-  tan 


!iii=  1  J* 

Z£h  2a3^ 


tan 


(10) 


where 

(i.  s  the  gravitational  constant  of  the  Moon 
fr  =  the  mass  magnitude  of  the  i**1  disk 
a  =  the  radius  of  the  disk 
By  solving  the  quadratic  equation  (Ref.  3) 


K2  +  (a2  -  R2)K  -  z2h  a2  =  0 


(ID 


for  K,  and  choosing  the  positive  root,  its  value  is 


a2)  + 


2.2  ^  .  2 
a  )  +  4  z 


ih 


(12) 


where 


R2  =  |[M)  r  -  P.|2 


and  where  £  is  the  selenocentric  position  vector  of  the  vehicle.  The  quantity 
[M]  is  the  transformation  matrix  from  the  mean  selenocentric  frame  to  the 
selenographic  frame. 
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Transformation  of  the  above  accelerations  to  the  selenographic  frame  is 
defined  as 


£i  =  R3  WR2(|-  *)h 


Thus,  acceleration  due  to  disk  masses  imbedded  in  the  Moon  becomes 


1 


The  differentially  correctable  disk  mass  parameter  is  the  mass  ratio  p.,  the 
variational  equation  for  which  has  the  nonhomogeneous  term 


where 


-o  rM-i]  — i 

r=  r  J  -3jrr 


3D.  T  rp  .  ,  3r 

•^T  =  R3  MrI(z-*)-5T 

l  r 1 


(aK1/2)/(a2+K)  -  tan"1  (aK"1/2)j  x£h 


-^  =  |pa“3  [(aKi/2)/(a2+K)  -  tan"1  (aK"1/2)|  y£] 


.-1  ,^-Vz 


-1/2, 


tan  A  (aK  '  )  -  (aK  *'“)  z 
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The  derivatives  of  the  disk  mass  acceleration  with  respect  to  the  vehicle 
position  £  and  velocity  £  are  given  by  the  matrices 


and 


3r0 

“5T=  f°l 


where  I  is  a  3  X  3  matrix. 
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GRAVITATIONAL  ACCELERATION  DUE  TO  THE  MOON 


Acceleration  due  to  the  Moon  is  computed  from  the  gradient  of  the  potential 
function 


U  '  r  1  +£  (“f)  £Pn‘  (sin  5)  (Cnm  cos  +  S, 
n=l  m=0 


sin  mi 


where 


)  <18> 


p  =  gravitational  constant  GM  of  the  Moon 

r,6,\  =  selenographic  distance,  latitude,  and  East 
longitude  of  the  vehicle,  respectively 

a^  =  mean  equatorial  radius  of  the  Moon 

P™  =  Legendre  associated  function  of  the  first  kind,  of 
'  degree  n  and  order  m 

C  _ ,S  =  numerical  coefficients. 

nm  nm 

The  gradient  of  the  potential  U  is  computed  in  the  radial,  east,  and  north 
directions  of  the  local  horizontal  coordinate  system.  Thus 


r,  =  M_1J[R]VU 


where  V  is  the  gradient  operator,  and  [m_1]  is  the  rotation  matrix  from  the 
selenographic  (MF  true  equator)  frame  to  the  seienocentric  (MCI 
mean  equator  and  equinox)  frame.  Thus 


[M]  =  [SP] 

=  [S][P] 
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where 


[  P]  =[  P(£q>  Z,  9)]  precession 
where 


=  right  ascension  of  the  mean  celestial  pole  of 
°  date,  referred  to  the  mean  equator  and  equinox 
at  the  base  date 

ir+Z  =  right  ascension  of  the  mean  celestial  pole  at  the 
base  date,  referred  to  the  mean  equator  and 
equinox  of  date 

0  =  tt/2  -  declination  (i.e.,  the  north  polar  distance) 
of  the  mean  celestial  pole  of  date,  referred  to 
the  mean  equator  at  the  base  date 

and 

[S]  =  (S(fl',A,i,  *,?,-)) 
where 

n',A,i  =  lunar  Eulerian  angles 
-  nutation  in  longitude 

7  -  true  obliquity  of  the  ecliptic  to  the  Earth's 
equator 

T  =  mean  obliquity  of  the  ecliptic  to  the  Earth's 
equator. 

Note  that  [  M]  (see  Ref.  6)  transforms  the  mean  equator  from  MCI  to  the  true 

MF  coordinates.  The  rotation  matrix  [  R]  in  Eq.  (19)  is  defined  as 


[R] 


cos  6  cos  \ 
cos  5  sin  \ 
sin  6 


-sin  X.  -sin  6  con  \ 
cos  \  -sin  6  sin  \ 
0  cos  6 


which  is  the  transformation  from  the  local  horizontal  frame  to  the  selenographic 
frame . 
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The  components  of  the  VU  are 


9U 

*1  =  TZ 


r  00  /a  \n  n 

+  !>(-?)  2C(s“6>*= „m  cos  mX  +  Snm  sin  mX) 
r  L  n=l  m=0 


-  1  9U 

“  r  cos  6  "5\ 


r  cos  6 


2(“r)  £niP™(sin6HC 


sin  rrA  -  S  cos  m\) 
nm  nm 


1  9U 

=7"36 


w  ,  v 

^  ^  \  r^J  ^  ]  P™  (sin  8)  cos  6  (Cnm  cos  mX.  +  Snm  sin  m\)  (21) 


where  P^1  (sin  6)  is  the  derivative  of  the  Legendre  function  with  respect  to 
sin  6,  and  subscripts  1,2,3  refer  to  radial,  east,  and  north  directions. 

The  recursive  formulas  used  in  computing  the  Legendre  associated  functions 
and  their  derivatives  are  divided  into  two  groups.  The  terms  with  m  =  0 
(zonal  harmonics)  are  computed  by  the  recursive  formulas 

P°sP  =  [  (2n  -  1)  sin  6  P  .-(n-l)P  J  /n 
n  n  L  '  n- 1  n-2J 


P'  cos  5  =  sin6  cos6  P;  .  +  n  cos6  P  , 
n  n- l  n- 1 


with  the  initial  values  Pq  =  Pj  =  1  and  Pj  =  sin6. 


The  terms  with  m  i  0  (tes serai  harmonics)  are  computed  by  the  recursive 
formulas 


T->rn  /  <■ 

P  /cos  6 
n 


=  j^(2n 


1)  3in  6  Pm  .  /cos  6  -  (n  +  m 
n- 1 


1)  Pm,/cos  5 
n -c 


/(n  -  m) 


and 


Pm  cos  6  =  (n  +  m)  Pm  .  /cos  6  -  n  sin  6  Pm/cos  6 
n  n- 1  n 

with  the  initial  values  Pm  ,  =  0,  and  Pm/cos  6  =  1*3  •  •  •  (2m- 1)  cos  m“*6. 

There  are  two  normalizations  possible  in  TRACE66  (Ref.  7)  when  the  Cnm 
and  Snm  coefficients  are  specified.  The  unnormalized  form  is  used  in  the 
above  potential  expression. 

Variational  equations  for  the  differentially  correctable  gravitational  param¬ 
eters  p,  Cnm,  and  Snm  are  obtained  by  differentiating  the  relations  in 
Eq.  (21).  These  derivatives  constitute  the  nonhomogeneous  terms  of  the  p, 
Cnm.  and  Snm  variational  equations 
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/  ,  f  co  /a  v“  n 

if  *7  *+£("  + *)(■?)  X>“ 

r  L  n=J  m=0 


(sin  6)(C  cos  mi.  +  S  sin  mX.) 


-^r  =  _ - - 

r  r  cos  6 


r  ^  n 

.  °o  /a  \  n 

-  /  (~~)  mP^1  (sin  5)(Cnm  sin  mX.  -  cos  mi) 

os  6  4  1 

_n=  1  m=0  -I 


CO  /  -*  .  li 

Ef'““ 


6)  cos  6  (C  cos  mX.  +  S  sin  m\) 
nm  nm 


nm  r  '  / 


in  6)  cos  mi 


=  -2-'e  (-f1)  mpr<sin6) 

r  cos  6 


sin  mX. 


a  n 

®3  {jl  /am\  -r^m'  ,  .  rl  c  N 

"9C —  =  2  ("T")  ^n  (sin  COS  ®  cos 
Lnro  r  L'  ' 


"3S  =  2 

nm  r 


—  =  (n  +  1 )  F™(sin  5 )  sin  ra'k 


d%2  v 

"53  =  2  H 

nm  r  cos  8 


in  5  )  cos  mi. 


[/a  'R 

—  =  — 1 f — ~*)  P™  (sin  6)  cos  6  sin  mi 
m  t  l* 


*35  "  2 

nm  r 
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Note  that  the  nonhomogeneous  terms  in  Eq.  (22)  have  the  general  form 


W 


8VU 


(23) 


where 


vu  =  [g1,g2,g3]T 

v  =  (J.,  C  ,  or  S 
1  n  nrn 


nm 


The  derivatives  of  with  respect  to  the  vehicle  position  £  and  velocity  r_  are 
the  matrices 


“ST 


[  R] 


du  9r 


=  (r]Ih][r]t 


(24) 


where  £  =  (g^,g2>g3)^  is  the  acceleration  in  the  local  horizontal  system,  and 
u  is  the  position  vector  in  the  horizontal  system  with  the  partial  derivative 
matrix  H  given  by 


-  g2 

9r 

cos  6 

~W  " 

*2 

%2/3\  +  g1  -  tan  6  g3 

*2 

9r 

cos  5 

~w 

%3 

9g2/9\ 

9g3, 

cos  5 

~w  ■ 

(25) 


and 


=  [o]. 
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3.3 


PLANETARY  GRAVITATIONAL  ACCELERATIONS 


Gravitational  acceleration  due  to  planetary  bodies  is  given  by 


(26) 


where 

=  mass  ratio  of  the  i**1  body  relative  to  the  Moon 
=  MCI  position  vector  of  the  i^1  perturbing  body 
£  =  MCI  position  vector  of  the  vehicle 

Note  that 

ECI  ECI 
r.  =  r.  -  r.  . 

—l  —l  — M 

where 

yECI  -  Earth-centered  inertial  (ECI)  position  vector  of  the 
ifck  perturbing  body 

£^<-'*  =  ECI  position  vector  of  the  Moon 

[Ephemerides  are  obtained  by  making  the  necessary  coordinate  frame  rotation 
and  interpolating  a  planetary  ephemeris  file  created  from  a  Jet  Propulsion 
Laboratory  n. aster  tape  (Ref,  8)], 

There  are  no  differentially  correctable  planetary  gravitational  parameters. 
With  respect  to  vehicle  position  £  and  velocity  £,  the  derivative  matrices  of 

h  are 


(r-r.)(r-r.)T 


where  I  =  the  3X3  identity  matrix,  and  (9£^)/(9£)  =  [o]. 


(27) 
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A 


3.4 


PARTIAL  DERIVATIVES 


If  Oc  denotes  the  vector  of  computed  measurements  and  P  the  parameter 
vector,  the  measurement  partial  derivatives  are 


90  90  /9r\  90  /  9r\ 

"5#  =  ^rl9p/+  ~dF  Isi/  (28) 


if  the  components  of  P  are  equation-of-motion-dependent  parameters.  Then 


are  called  geometric  measurement  partials 


are  called  trajectory  partials,  which  are 
solutions  of  the  variational  equations 


Note  that  for  measurements  obtained  from  Earth-based  sensors,  TRACE66 
considers 


9m^  9m^ 


9p.  “  "c5F 


a-ECI 

-ft - 

ECI  -MCI 


^XMCI  \ 

r^r) 


K  8^ECl/8-Md  \ 


*9r 


-ECI  8-MCI 


(29) 
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where 


3.4.1 


fegCI  _  5(-MCI  +  -M) 
9-mci  8-mci 


8-ECI 

S-MCI 


8^— MCI  +  — 
a-MCI 


=  I 


I  =  identity  matrix 


£^  =  ECI  position  vector  of  the  Moon 
m£  =  computed  jfch  measurement 

=  ifch  equation- of-motion  parameter 
Trajectory  Partiala 


In  TRACE66,  trajectory  partials  for  equation-of-motion-dependent  parameters 
are  obtained  by  numerically  integrating  the  vector  variational  equation 


8r  8i^  8r ,  Sr., 

■37  =  IF  +  “ST  +  “ST 


dr_  9r  ( £ ,  £,  y  ) 

+  - g- - 


(30) 


where  y  is  some  parameter  of  the  parameter  vector  £.  Thus 

9r 

The  equation-of-motion  parameters  available  in  the  Lunar  Analyzer  are 

Vehicle  State  *  -  {(£.£)0,  tQ} 

Point  Mass  - 

Disk  Mass  -  {p^} 

Central  Body  -  {p,  C  ,  S  } 

7  nm  nm ; 

^Initial  vehicle  state  parameters  may  be  expressed  in  various  coordinate 
systems  (Ref.  3). 
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3.4.2  Geometric  Measurement  Partials 

The  observation  measurements  used  in  the  determination  of  the  lunar  gravity 
field  consist  exclusively  of  JPL./NASA  Deep  Space  Network  long-count  two- 
way  doppler  measurements,  denoted  by  the  symbol  CC3.  Section  4  presents 
a  detailed  description  of  the  measurement  model  and  the  equations  used  to 
compute  the  geometric  measurement  partials,  8CC3/9i:and  3CC3/3K  Here 
it  suffices  to  say  that  these  partial  derivatives  are  independent  of  the  number 
and  nature  of  the  parameters.  They  are  computed  based  only  on  the  position 
and  velocity  of  the  satellite  and  the  earth-based  sensors  at  the  appropriate 
times . 
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4.  MEASUREMENT  MODEL 


Observation  measurements  used  in  the  determination  of  the  lunar  gravity 
field  are  the  JPL/NASA  Deep  Space  Network  long -count  two-way  doppler 
(see  Ref.  9).  A  stable  oscillator  (see  Figure  4)  is  used  to  generate  a  reference 
signal  that  is  transmitted  to  the  spacecraft.  The  reference  signal  frequency 
is  subtracted  from  the  return  signal  frequency,  which  results  in  the  doppler 
tone.  The  two-way  doppler  data  constitute  an  accumulative  counter  reading 
value;  the  integrated  cycle  count  divided  by  its  time  tag  difference,  in  seconds, 
becomes  the  two-way  doppler  observable.  The  corresponding  time  tag  for 
the  observable  is  the  time  halfway  between  the  counter  reading  time  tags. 

When  one  performs  the  computations  described  in  this  section,  it  is  necessary 
to  calculate  position  and  velocity  of  the  spacecraft  relative  to  the  earth  at 
dmes  that  may  not  correspond  exactly  to  an  integration  time.  Therefore,  it 
is  necessary  to  interpolate  as  well  as  to  translate  in  order  to  obtain  the 
desired  results.  The  Appendix  contains  some  notes  on  the  importance  of  the 
order  and  the  accuracy  of  these  computations. 
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32 


Figure  4.  Simplified  Two-Way  Doppler  System 


4.  i  LONG-COUNT  TWO-WAY  DOPPLER  FREQUENCY 

The  computed  two-way  doppler  observable  (see  Ref.  9)  is  expressed  as 


CC3  =  -rr-  (T 04. 25339367)  Ap  (31) 

ctC 

where 

Lj,  =  transmitter  frequency  —  22,  045,  600  Hz 
c  =  velocity  of  light  =*  299792. 5  km/sec 
t^  =  count  time,  usually  60  sec 

Ap  =  topocentric  range  difference 


Geometry  for  the  long-count  two-way  doppler  frequency  is  shown  in  Figure  5. 

4.2  COMPUTATIONAL  ALGORITHM  FOR  LONG- COUNT 

DOPPLER  FREQUENCY 

The  notation  used  in  the  formulation  for  long-count  doppler  frequency  is 

tg  =  time  at  which  frequency  count  begins  at  the  receiver 

tp  =  time  at  which  frequency  count  ends  at  the  receiver 

— B’— E  =  8eocen*:r^c  position  vector  of  receiver  R  at  time 
tg  and  tg,  respectively;  that  is,  Rg  =  R(tg), 

R„  =  R(t„),  where  R  is  the  position  vector  of  the  receiver. 

—Jc Li  —  Sit 


Ig  and  tp,  respectively,  minus  signal  transit  time  from 
spacecraft  to  receiver 

geocentric  position  vector  of  spacecraft  at  time  t  and  tOT,, 

SB  SE 

respectively 
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geocentric  position  vector  of  transmitter  T  associated  with 
beginning  and  end  of  count,  respectively 

total  count  time  =  t  -  f 

E  B 

time  tag  of  data  on  JPL  B2  tapes  (see  Ref.  10) 


a>e  =  Earth's  angular  rate 


P 


(  ) 


topocentric  range. 


The  computational  algorithm  consists  of  the  following  steps. 
Step  1 .  tg  and  t^,  are  obtained  from  t  and  t^  by 


(32) 


Step  2.  From  tfi  and  tg,  R_B  and  RE  may  be  found  directly 
by  known  formulas . 

Step  3.  £b  is  obtained  by  an  iterative  process  which  takes  into  account 
signal-transit-time  delay. 

a.  It  is  assumed  that  the  geocentric  position  ,r  of  the 
spacecraft  has  been  integrated  from  the  equations 
of  motion  and  is  a  known  function  of  time  that  may 
be  interpolated  to  the  required  accuracy  at  any  given 
epoch. 
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b.  First  estimates  of  the  transit  time  Atg^g  along  the 
SRB  leg  (Figure  5)  and  the  time  t^g  are  given  by 


At 


|£bVI 


SRB 


1 1  =  t  -At*' 

CSB  lB  SRB 


(33) 


c.  The  initial  estimate  of  the  range  is 

PSRB  =  l-B^SB^  "  ~ ttB)  I  +  ApSRB 


(34) 


where  is  a  refraction  correction  calculated 

once  and  held  fixed  for  succeeding  iterations, 

d.  An  improved  estimate  of  transit  time  is 


At 


PSPJ 


SRB 


(35) 


from  which  is  obtained 


2  2 
fcSB  =  fcB  "  AtSRB 


(36) 


e.  Step  3.c  is  repeated  to  find  that 


2  2 

PSRB  =  I-B^Sb)  "  — (fcB)  f  +  ApSRB 


(37) 
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f.  Convergence  is  obtained  if 


Ipsrb^sb)  "  Psrb^sb)  | 


<  10~7km 


Step  4 . 


Otherwise  return  to  Step  3.d. 

After  tjjig  and  R.(t^gj  are  found,  _Tg  is  determined.  Thus,  by 

a  procedure  similar  to  the  above  and  by  starting  with 

t  -  tn  Atn 
CTB  ZSB  'SRB' 


(38) 


where  both  quantities  on  the  right  side  are  obtained  from 
the  last  iteration  of  Step  3.  The  first  estimate  of  Pg^g  is  then 


PSTB  =  lr(fcSB)  '  — B^Tb)!  +  ApSTB 


At 


PSTB 


STB 


(39) 


and  ApgTB 


is  the  refraction  correction 


*TB  =  fcSB  “  AtSTB 

PSTB  =  I-^Sb)  "  —B^TB)  1  +  ApSTB  (40) 


Step  5. 


Steps  3  and  4  are  repeated  to  solve  for  PgRE  and  PS^E*  For 

successive  observations  separated  by  a  time  interval  t~, 

o 

th 

PSRE  an<*  PSTE  m  n  observation  are  same  as 
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PSRB  and 

For  observations  separated  by  more  than  t^,,  the  calculations 
of  p  must  be  made  for  all  four  legs. 

Step  6.  The  refraction  correction  for  each  leg  is  computed  from 

.  .0018958N  ,  .  „  .  ,n^-1.4 

Ap  = - “) -  (sm  E  +  •  06483) 


pgTB»  respectively,  of  the  (m  +  l)s*-  observation. 


Step  7 . 


where  E  is  the  elevation  angle  of  the  receiver  (or 
transmitter)  on  the  appropriate  leg  at  tgg  (or  t,pg,  etc., 


depending  on  the  leg);  and  N  is  the  surface  refractivity 

in  ui.its  of  10'6,  nominally  equal  to  300.  The  correction 
is  added  to  each  calculated  p^  y 


Form  the  quantity  Ap  (including  refraction  corrections)  as 


A  p  =  PSTE  +  PSRE  -  pST3  -  pSRB 


Step  8 .  Calculate  a  two-way  doppler  from 

fT 

CC3  = -~  (104.25339367)  Ap  (41) 

ccC 

4.3  MEASUREMENT  PARTIAL  DERIVATIVES 

The  measurement  partial  derivatives  9CC3/3p  are  calculated  from  the  chain 
rule  and  simplifying  approximations 


PSTE  '  PSTB  =*  pSTEfcC 
PSRE  "  PSRB  =<  pSREfcC 
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where 


(p  P)STE  =  (-E  "  -E),(-E  "  ^E* 

P^SRE  =  E  "  — E^ *  — E  "  — E^ 
Thus,  Eq.  (41)  becomes 


CC3  =  B[  pSTE  +  pSRE] 


(42) 


where 

fT 

B  =  —  (104.25339367) 

Let  the  components  of  the  vector  £,  R,  TT  be  expressed  as 


T 

r  =  (x, y,  z) 

T 

~~  ZR^ 

T  ^  (xT,yT,  zt)T 


then 


3CC3  _ 
~^T  =  B 


17*  +  “eyT  -  P  Lx\  +  (; 

V  P  'STE 


+  u  yn 
e7  R 


(43) 


In  Eq.  (43),  the  expression  (  )gj£  means  that  all  quantities  inside  the  paren¬ 
theses  are  evaluated  along  the  leg  STE;  and  L  is  a  direction  cosine  defined  in 


(L)STE  =  (LV  Ly'  LzWe 


TE  =  ( 


x  -  x„ 


y  -  Y' 


z  -  z 


A 


STE 


(44) 
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Also, 


9x  “  B[(Lx)STE  +  ^-x^SRe] 

=  B[(Ly)STE  +  (Ly}SRE] 

3z  =  B[(Ljz)STE  +  (Lz  *SRe]  (45) 


Now  let  p  be  any  parameter  to  be  estimated,  and  S  =  (r,jr)  where  r  =  (x,y.  z)T 
T 

and  £  =  (x,y,z)  .  Then  the  6X1  matrix  of  solutions  to  the  variational  equa¬ 
tions  3S/ 3p  -  9(x,  y,  z,x,  y,  z )/ 3p  is  obtained  by  integrating  a  set  of  differential 
equations,  as  described  in  Section  3.  Therefore,  the  required  measurement 
partial  derivative  needed  for  the  differential  correction  is 


3CC3  _  / 3CC3 \  /9-\ 
3p  "  \  3S  j  \~5pj 


(46) 
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5.  NORMAL  MATRIX  ACCUMULATION 


5. 1  NORMAL  MATRIX  ACCUMULATION  PER  ARC 

The  augmented  normal  matrix  is  a  weighted  matrix  product  of  the  form  A^WA, 
where  A  is  the  matrix  of  measurement  partial  derivatives  and  W  is  the  inverse 
of  the  measurement  error  variance- covariance  matrix.  The  jth  rQ\v  of  the  A 
matrix  has  the  form 

9m^  9m^ 

_ £  _ c_  _ c  nj 

9P,  »  DpT  •  *  *  *  »  ~W  *  me 

12  n 

where 

mJ,  =  the  computed  measurement 
Pj,  =  the  i^h.  parameter,  i  =  1,  .  .  .  ,  n 

=  the  partial  derivative  of  the  j  1  measurement  with 
1  respect  to  the  ifck  parameter 

OJ  =  the  jtk  measurement  residual;  O^  =  m3  -  m3  , 
me  J  .  me  c 

where  m.J  is  the  actual  measurement 

Let  Oc  denote  the  vector  of  computed  measurements.  Then,  when  the  param¬ 
eters  are  quantities  appearing  in  the  equations  of  motion  (including  initial  con¬ 
ditions),  the  measurement  partial  derivatives  are  computed  from  the  chain 
rule 


90  90  3r  90  af. 

c  c  —  ,  c  Q£. 

DF  =  ~5r"  5p  +  ~W~  c5p  (47) 

where  9r/ 9P  and  9ry  9P  are  referred  to  as  trajectory  partials,  and  are  solu¬ 
tions  of  the  variational  equations.  The  matrices  90c/9r;  and  90c/3£,  and  the 
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columns  of  A  that  correspond  to  parameters  not  appearing  in  the  equations 
of  motion  (station  locations  and  measurement  error  terms)  are  referred  to 
as  the  geometric  partials  and  are  computed  analytically. 

If  W  were  a  full  matrix,  it  would  be  necessary  to  compute  the  entire  A  matrix 
at  one  time  in  order  to  perform  the  matrix  multiplications.  However,  in  the 
Lunar  Gravity  Analysis,  it  is  assumed  that  all  observations  are  independent 
and  have  the  same  variance.  Therefore,  W  is  of  the  form 

W  =  1  /cr2  I 

where  I  is  an  identity  matrix.  The  normal  matrix  may  be  written  as 

ATWA  =  l/o-2  ATA 

=  ‘/"2  £>„»*]  <«> 

i=l 


where 


a..  =  the  i 
ij 


;th 


element  of  the  matrix  A 


^2i>- 


i=  1 


a^  =  the  j  -  k**1  element  of  the  matrix  A^WA 


The  first  n  columns  of  the  A^WA  comprise  the  normal  matrix  of  the  least- 
squares  problem.  The  last  column  is  the  right-hand  side  of  the  system 
A^WO  with  last  element  of  that  column  equal  to  the  weighted  sum  of  the 
squares  of  the  residuals  0^cW  Omc«  Note  that  by  using  Eq.  (48), 
the  A  matrix  needs  to  be  formed  only  one  row  at  a  time,  and  the  contribution 
from  that  particular  observation  is  summed  into  the  A  WA.  Furthermore, 
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since  the  normal  matrix  is  symmetric,  the  accumulation  is  done  only  for 
k  >  j,  thereby  eliminating  the  computations  for  all  elements  below  the 
principal  diagonal. 

5.2  MERGING  THE  SINGLE-ARC  MATRICES 

A  computer  program  external  to  TRACE66  takes  the  single-arc  normal  matrices 
generated  as  described  above  and  merges  them  into  one  large  linear  system. 

All  input  and  output  matrices  of  the  merge  procedure  are  stored  on  magnetic 
tape.  A  priori  inverse  variances  may  be  added  to  the  diagonal  elements 
corresponding  to  mass  parameter  magnitudes  at  the  time  of  the  merge. 

Vehicle  state  parameter  a  prioris  can  be  added  if  desired  at  the  time  of  the 
single-arc  matrix  accumulation  in  TRACE66.  The  form  of  the  final  merged 
normal  matrix  is  shown  in  Figure  6. 
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MASS  MAGNITUDE 
PARAMETERS 


NONZERO  ELEMENTS  CORRESPOND  TO 
MASS  PARAMETERS  SELECTED  FOR 
CORRESPONDING,  VEHICLE 
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Figure  F.,  Form  of  the  Merged  Normal  Matrix 


6,  SOLUTION  OF  THE  NORMAL  EQUATIONS 


Two  alternate  techniques  are  used  to  solve  the  normal  equations.  The  first 
method,  for  problems  with  a  limited  parameter  set  (i.e.  ,  a  maximum  of 
100  parameters),  computes  the  inverse  of  the  normal  matrix  as  well  as 
the  parameter  corrections.  In  addition,  a  side  condition  can  be  imposed 
that  bounds  the  magnitude  of  the  parameter  corrections.  This  capability, 
including  the  solution,  is  completely  self  -contained  in  TRACE66. 

The  second  method  is  used  for  problems  v/ith  an  effectively  unlimited 
number  of  parameters.  It  uses  mass  storage  devices  during  the  solution 
process,  and  therefore  does  not  require  the  retention  of  the  entire  normal 
matrix  within  core.  With  this  method  the  solution  of  the  normal  equations 
is  performed  in  a  program  exterior  to  TRACE66.  (Note  that  the  inverse 
of  tne  normal  matrix  is  not  computed  in  this  method.) 
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SOLUTION  .FOR  A  LIMITED  PARAMETER  SET 


6.  1 

The  complete  solution  procedure,  including  the  bounding  feature,  is  fully 
discussed  in  Reference  11,  and  will  be  only  briefly  described  here.  The 
reference  states  that  if  a  function  to  minimize  is  given,  as  well  as  a 
bounding  condition  of  the  form 

f  (AP)  =  ^AAP  -  OmcjTW(AAP  -  Omc)  (49) 

and 

(GAP)T  (GAP)  <  1 

the  value  of  AP,  which  minimizes  f  subject  to  the  constraint,  is  the  solution 
of  the  linear  system 


(ATWA  +  ZGTG)  AP  =  ATW  O  (50) 

me  v  ' 

where  Z  is  a  Lagrange  multiplier.  The  proper  value  of  Z  may  be  deter¬ 
mined  by  an  iterative  process.  Therefore,  the  problem  is  always  reduced 
to  the  finding  of  the  solution  to  a  system  of  the  form 


CX  =  B 


where  C  is  an  n  x  n  positive  definite  symmetric  matrix:  X  is  an  m  vector 
of  unknowns;  and  B  is  an  m  vector  referred  to  as  the  right-hand  side.  The 
procedure  used  to  compute  X  and  invert  the  matrix  C  is  finite  (noniterative) 
and  applicable  only  to  symmetric  matrices. 

It  is  based  upon  the  fact  that  a  symmetric  matrix  can  be  decomposed  into 

T 

a  product  of  the  form;  C  =  LDL  where  L  is  lower -triangular  with  -i  on 
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the  diagonal;  and  D  is  a  diagonal  matrix.  In  such  a  representation, 
det(L)  =  ±  1,  and  det(D)  =  det(C).  Therefore,  L  *  exists,  and  is  also 
lower -triangular  with  -1  on  the  diagonal;  and  D  has  no  zero  diagonal  element 
if  C  is  non-singular.  Thus 

D  =  L_1  C  (LT)-1  =  SCST 

C-1  =  (LV1  D"1  L-1  =  ST  D"1  S  (51) 

where  S  =  L  *. 

The  procedure  builds  S  one  row  at  a  time  with  a  bordering  technique. 

Since  D  is  diagonal,  once  S  is  obtained,  the  calculations  of  X  and  C  *  are 
reduced  to  matrix  multiples. 

6.  2  SOLUTION  FOR  AN  UNLIMITED  PARAMETER  SET 

When  a  great  many  unknowns  are  being  considered,  the  system  is  solved 
by  using  a  technique  based  upon  the  Gaussian  elimination  method.  The 
procedure  is  to  successively  eliminate  variables  from  the  system  until  the 
problem  is  reduced  to  one  equation  with  one  unknown.  By  using  the  solution 
of  this  equation  as  a  starting  point,  values  for  the  remaining  unknowns  are 
obtained  through  back  substitution. 

Symbolically,  the  process  is  expressed  as  follows:  Let  an  n  X  n  system 

be  of  the  form  CX  =  B,  where  C  is  ?n  n  X  n  matrix  [c . . ] ,  X  an  m  vector 

of  unknowns  [x.J,  and  B  an  in  vector  [b^].  Then  eliminate  the  k  column 

and  row  (k  =  1,  n  -  1)  by  modifying  all  c..  and  b.,  where  i,  j  >  k  in 

ij  i 

the  manner 
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(52) 


c(.  =  C..  - 
1J  ^  Ckk 


c.,  b, 

b(  =  b.  -  ^  k 

1  1  Ckk 


Finally,  after  (n  -  1)  eliminations,  the  system  will  be  reduced  to 


c  x  =  b 
nn  n  n 


where  c  and  b  have  been  modified  (m  -  i)  times.  Then, 
nn  n  '  '  ’ 


x  =  b  /c 
n  n  nn 


Using  the  solutions  for  x^+j>  •••  ,  x  ,  produces  the  solution  for  x^,  which 
is 


where  the  c's  and  b's  are  the  modified  coefficients  after  all  eliminations. 

A  slight  modification  of  the  technique  can  be  used  to  take  advantage  of  a 

symmetric  matrix.  As  the  symmetry  property  can  be  expressed  as 

c..  =  c..  for  all  i  and  j,  the  basic  modification  eciuation  can  be  rewritten  as 

1  i  vi  J 


c'. . 
U 


c.. 

U 


Ckk 


(54) 
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The  above  needs  to  be  calculated  only  for  all  j  >  k,  and  for  all  i  >  j  for 
each  j.  This  eliminates  the  computation  of  all  modified  elements  above  the 
principal  diagonal  of  the  matrix. 

Another  sophistication  of  the  technique  significantly  increases  the  algorithm's 

efficiency.  It  reduces  the  number  of  multiples  required  when  the  matrix 

is  sparse  with  strings  of  consecutive  zeros  in  some  columns.  When  c., 

til 

is  zero,  none  of  the  elements  in  the  i  row  need  be  modified  for  eliminating 
th. 

the  k  unknown.  By  using  an  indexing  scheme  that  points  to  the  beginning 
and  end  of  strings  of  zero  elements  in  the  matrix,  computation  time  is 
considerably  reduced. 

For  systems  of  the  size  being  considered,  a  large  number  of  central 
memory  words  of  data  storage  would  be  required  to  obtain  the  solution  in 
core.  For  example,  a  system  on  the  order  of  2000  would  require  in  excess 
of  2,  000,  000  words.  Because  this  amount  of  core  storage  is  not  feasible, 
a  technique  of  file  usage  (either  magnetic  tape  or  disk)  was  devised  to  take 
advantage  of  mass  storage  devices.  Three  files  are  used  with  the  following 
designations  and  purposes: 

ATAPE  Contains  the  matrix  C  stored  as  lower 

triangular  by  columns 

TTAPE  Contains  the  transformations  necessary 

for  eliminating  each  column,  along  with 
start/stop  pointers  for  strings  of  zero 
elements  in  that  column 

NBACK  Contains  the  transformed  A  matrix  to  be 
used  for  back  substitution. 

A  buffer  of  the  maximum  available  size  is  provided  for  the  program.  The 
number  of  file  manipulations  required  is  inversely  proportional  to  the  size 
of  this  buffer.  Therefore,  an  increase  in  buffer  size  results  in  a  reduction 
of  job  throughput  time.  The  right-hand  side  (B  vector)  is  kept  in  core. 

The  buffer  is  filled  from  ATAPE  with  as  many  columns  of  C  as  will  fit. 
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Starting  with  the  first  column  of  C,  which  i3  in  the  buffer,  and  working 
sequentially  towards  the  last  column  in  it,  transformations  are  computed; 
applied  to  the  remaining  columns  in  the  buffer  and  to  the  B  vector;  and  then 
written  on  TTAPE.  When  the  last  column  has  been  prDccssed,  the  contents 
of  the  buffer  are  written  on  NBACK.  The  buffer  is  refilled,  TTAPE  is 
rewound,  and  the  transformations  are  read  in  and  applied  to  all  columns. 
The  process  is  then  repeated  for  this  buffer  load,  starting  with  the  com¬ 
putation  of  the  transformation  for  the  first  column.  (Note  that  additional 
transformations  are  added  to  TTAPE  for  every  buffer  load.  )  After  the 
last  column  of  the  matrix  has  been  processed,  back  substitution  is  begun. 
The  necessary  transformed  C  matrix  is  obtained  by  initially  reading  the 
last  buffer  load  from  NBACK,  processing  it,  backspacing  NBACK  twice, 
reading  the  next-to-last  buffer  load,  and  so  on,  until  the  full  solution  is 
obtained. 

The  total  execution  time  depends  upon  the  order  of  the  system  being  solved, 
the  size  of  the  buffer  supplied,  and  the  type  of  mass  storage  devices  being 
used.  Central  processor  (CP)  time,  however,  is  only  slightly  affected  by 
buffer  size  or  storage  device.  The  following  empirical  equation  (although 
slightly  pessimistic  for  large  systems,  n  >  1000),  has  proved  reliable 
for  typical  problems: 


CP  0.  2  +  0.  0075  n  +  .  00003  n2  +  1.  77  X  10'7  n3  (sec)  (55) 

As  successive  columns  of  the  matrix  are  eliminated,  divisions  are  performed 

by  the  current  diagonal  element  of  that  column  (i.e. ,  the  one  transformed 

by  all  previous  eliminations).  Because  the  matrix  is  symmetric  and  positive 

definite,  all  transformed  diagonals  are  smaller  than  the  original  diagonal 

element.  A  comparison  of  the  element  used  for  the  division  with  the  original 

diagonal  element  provides  information  regarding  the  conditioning  of  the 

matrix.  If  the  ratio  of  these  numbers  is  less  than  some  i  (for  example, 

-5 

£  -  10  ),  then  the  corresponding  unknown  is  said  to  be  pooily  determined. 
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On  option,  that  unknown  may  be,  in  effect,  eliminated  from  the  system 
by  zeroing  out  the  entire  row  and  column,  and  setting  the  diagonal  element 
to  a  large  number  (for  example,  Such  a  procedure  enhances  the 

stability  of  the  system.  In  the  case  of  the  Lunar  Gravity  Field  determina¬ 
tion,  the  solution  is  essentially  unaltered  for  a  small  number  of  parameter 
eliminations. 
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7.  COMPOSITE  LUNAR  GRAVITY  FIELD 


Several  representations  of  the  lunar  gravity  field  have  been  obtained  recently 
(see  Ref,  12).  The  fields  differ  more  in  the  methods  of  analysis  than  in  the 
basic  data,  since  all  the  results  were  based  on  the  doppler  tracking  of 
Lunar  Orbiters.  A  procedure  for  combining  a  rnascon  gravity  model  (see 
R  13)  with  spherical  harmonics  to  obtain  a  composite  model  of  the  lunar 
gravity  field  is  described  in  Sections  8.  1  and  8.  2. 

7.  1  SYNTHESIS  OF  THE  COMPOSITE  FIELD 


Let  the  surface  of  a  sphere  of  radius  r  be  divided  into  two  regions  (see 
Figure  7).  Region  R^  is  the  one  under  which  mascons  have  been  obtained, 
and  R2  is  the  remainder  of  the  sphere.  Now  let  a  network,  or  grid  of  points 
on  the  sphere,  have  coordinates  (r,  <Jy  A^)  at  the  jth  point.  Let  be  the 
number  of  points  in  Rj,  and  j2  the  number  in  R^;  and  J  the  total  number, 
or  J  =  jj  +  j2«  At  each  of  the  points  in  R^,  the  vertical  component  of  the 
gravity  disturbance  gf  is  calculated  from  the  rnascon  distribution,  while 
in  R2,  gr  is  calculated  from  a  spherical  harmonics  model.  The  value  of 
r  is  chosen  so  that  the  altitude  is  100  km  above  the  mean  lunar  surface, 
because  most  of  the  best  doppler  tracking  for  the  rnascon  model  occurred 
at  this  altitude. 


1 1"1 

The  disturbance  potential  at  the  5  grid  point  caused  by  the  i  rnascon 
(see  Figure  8)  is  given  by 


V.. 


,Ai  f  2  2  ?  I' 1/2 

— —  =  p.  a  +  r  -  2ar  cos  vi;.. 

P- •  i 


(56) 


Preceding  page  blank 
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whore 


cos  Jj..  =  sin  d.  sin  $.  f  cos  <b.  cos  $.  cos  (A.  -  X.) 

Yil  i  J  i  J  J  i 

[x.  =  magnitude  of  i**1  rnascon 

a,  ,  X.  s  mean  lunar  radius,  selenographic  latitude, 
1  1  and  longitude,  respectively,  of  i*h  rnascon 

r,  A.  :=  selenographic  distance,  latitude,  and 
longitude  of  grid  point. 


The  corresponding  radial  component  of  gravity  disturbance  is 


9V. . 

_ 11 

9r 


“i  , 

“  T  (r 

pij 


a  cos  dj. .) 


(57) 


If  I  is  the  total  number  of  mascons  in  the  distribution  and  the  net  value  of  g„ 

r 

is  a  superposition  of  all  the  contributions,  then 


gT  (* 


•  V  V  ■  I  (•  5f)  ■ 


I 

E 


(r 


a  cos  ip. .) 


(58) 


ij 


Equation  (58)  is  then  used  to  obtain  a  gr  value  for  each  of  the  grid  points 
in  R^. 


At  each  of  the  j ^  points  in  I^,  a  gf  value  can  be  calculated,  from  the  radial 
derivative  of  the  disturbance  spherical  harmonics  work  function  (negative 
of  potential)  of  degree  N,  as 
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Figure  8.  Geometry  of  Grid  Point 
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N  n 


V(r,  4,  A)  = 


n=0  m=0 


Pnm  (9in  ,Cnm  008  mA  +  Snm  sin  mA) 


where 


/  \H+1 

=  g  a  y  V  (C  Yc  +  S  Y£ 

6o  \  a  I  nm  nm  nm  r 


5  lunar  mass 


r,  $,  A  s  spherical  coordinates  of  grid  point  in 

Hnm,  ^nm  =  normalized  coefficients 

P  =  associated  Lengendre  polynomial  (fully 
nm  normalized) 

So  51*m/a2 


The  radial  component  of  gravity  disturbance  being 


gr(r,  A)  =  -g02  2  (t) 

n  m  \  / 


(n  +  1)  (C  YC  +  S  Y9  )  (60) 

nm  nm  nm  nm7  ' 


a  gf  value  for  each  point  in  can  then  be  calculated  from  Eq.  (60). 

7.2  HARMONIC  ANALYSIS 

A  harmonic  analysis  of  the  radial  gravity  field  can  now  be  obtained  from 
Eqs.  (58)  and  (60).  Let  the  expansion  of  the  composite  field  g^  have  the  f  m 
of  Eq.  (60);  that  is 


f 


(61) 


gr(r»  A.) 


-  Sr 


EE(C 


(n  +  1)  (C  YC  +  S  YS  ) 
nm  nm  nm  nm 


n  m 


where  the  symbols  have  similar  meanings  and  p  is  the  degree  of  the  new 
expansion.  If  all  of  the  grid  values  of  gf  (denoted  g°)  are  assigned  equal 
weights,  the  coefficients  in  the  expansion  of  Eq.  (61)  are  easily  obtained 
from  the  orthogonality  properties  of  Legendre  polynomials,  so  that 


7.3 


AUXILIARY  CALCU LATIONS 


After  the  harmonic  analysis,  a  number  of  auxiliary  calculations  are 
obtainable,  such  as  the  following. 

a.  An  RMS  fit  over  regions  Rj  and  is  obtained  by 


r 


-.1/2 


RMSi=l2^ 


-i/2 

h 


rms2  = 


j-j^  -jl/2 

Z  Wj  (gr  - 

U=1 


.—1  /2 


(63) 


,1/2 


r  2  21  *  ** 

RMS  =  RMS^  +  RMS2  the  RMS  power  in  the  field 


-57- 


b. 


For  r  =  a,  the  power  spectrum,  or  degree  variance, 
of  the  potential  coefficients  is 


KMSt 


(64) 


c.  The  spectrum  of  vertical  gravity  at  r  =  a  +  h,  where 
h  =  100  km,  is  given  by 


2 

or 


n 


(65) 


d.  The  potential  coefficients  Cnm,  Snm  and  their 
unnormalized  equivalents  are  obtained  by 


(66) 


7. 4  REMARKS 

The  procedure  described  above  for  obtaining  a  composite  lunar  gravity  field 
is  similar  to  the  one  applied  to  the  harmonic  analysis  of  terrestrial  gra¬ 
vimetry.  In  the  latter,  measured  surface  gravity  is  assigned  to  squares 
with  measurements,  and  gravity  from  a  satellite  field  is  assigned  to  squares 
without  measurements.  However,  the  basic  data  base  is  completely  different 
for  the  two  cases,  as  the  expansion  obtained  at  100  km  can  be  analytically 
continued  to  the  lunar  surface  since  there  are  no  intervening  masses. 
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The  maximum  number  for  J  should  be  about  1640  corresponding  to 
5X5  deg  squares,  and  the  order  of  the  expression  should  have  a  limit  of 
20,  which  implies  the  solution  of  a  system  of  order  441. 

The  composite  lunar  gravity  field  computations  are  performed  outside  of 
TRACE66. 


8.  LUNAR  GRAVITY  ANALYSIS  FEATURES 


In  addition  to  the  normal  matrix  accumulation  and  least-squares  simulation 
options  described,  the  following  features  are  available  in  TRACE66. 

8.  1  MASS  PARAMETER  SELECTION 

The  mass  parameters  to  be  estimated  are  selected  by  the  following  criterion: 
The  acceleration  arising  from  a  given  mass  (point  or  disk)  integrated  over 
the  counting  interval  must  exceed  the  noise  threshold  of  0.  01  Hz  at  some 
point  along  a  given  trajectory  before  that  mass  is  selected  for  regression 
from  the  doppler  data  on  the  trajectory.  The  criterion  can  be  translated 
into  an  approximate  condition  on  the  distance  P  from  the  center  of  the 
point  or  disk  mass  to  the  satellite  trajectory;  that  is,  P  must  be  less  than 
some  minimum  distance. 

The  computational  algorithm  for  mass  parameter  selection  is 


a.  P=i[S]Tr-P.! 

b.  P  <  (»  +  ph)  (67) 

where 

P.  =  selenographic  position  vector  of  the  i**1  mass 

[s]  =  transformation  matrix  from  selenographic 
coordinates  to  selenocentric  frame 

a,P  =  altitude  bias  and  scale  factors,  respectively 

r_  =  selenocentric  position  vector  of  the  vehicle 

h  =  vehicle  altitude. 


Preceding  page  blank 
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8.2 


RESIDUAL  CALCULATIONS 


Calculated  and  plotted  at  each  observation  time  are:  the  filtering  of 
observation  data  blunder  points;  the  doppler  residuals;  spacecraft 
selenographic  latitudes,  longitudes,  and  heights. 
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APPENDIX 


INTERPOLATION  ACCURACY  IN  TRACE66  WITH 
APPLICATION  TO  LUNAR  EPHEMERIDES 

A.  1  METHOD 

TRACE66  uses  a  Hermitian  interpolation  scheme  to  calculate  positions  and 
velocities  at  times  not  exactly  equal  to  the  time  of  some  integration  point. 

The  Hermitian  polynomial  used  to  calculate  position  is  a  function  of  position, 
velocity,  and  acceleration  at  two  different  times,  one  on  each  side  of  the 
desired  unknown  time  point.  These  six  pieces  of  information  allow  the  deriva¬ 
tion  of  a  fifth-order  interpolating  polynomial.  The  formula  for  the  velocity 
is  the  derivative  of  the  position  equation,  and  is  therefore  a  fourth-order 
polynomial. 

A.  2  ACCURACY  OF  THE  METHOD 

The  accuracy  of  these  polynomials  depends  on  three  things: 

•  The  time  interval  or  step  size  between  the  known 
points 

•  The  relative  position  of  the  unknown  point  in  this 
interval 

•  The  nature  of  the  function  or  dependent  variable  to 
be  calculated;  i.e.,  sizes  of  the  function's 
derivatives 

Note  that  calculation  of  a  dependent  variable  in  the  exact  middle  of  an  interval 
leads  to  the  largest  interpolation  error  that  will  occur  in  that  interval.  All 
accuracy  figures  given  in  this  Appendix  correspond  to  interpolation  at  the 
interval's  midpoint  for  the  ephemeris  of  a  typical  ’ow-altitude  earth  satellite. 
The  results  of  an  empirical  study  of  interpolation  accuracy  versus  step  size 
are  presented  in  Fig.  A-i.  The  fifth-order  position  and  fourth-order  velocity 
curves  were  calculated  from  the  equations  used  in  TRACE66.  The  third- 
order  Hermitian  curves  are  plotted  for  comparison  purposes.  It  is  expected 
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that  a  smoother  function  such  as  the  ephemeris  of  a  high-altitude  earth  satel¬ 
lite  would  yield  greater  accuracy,  whereas  a  more  radical  function  might 
yield  somewhat  poorer  accuracy. 

A.  3  APPLICATION  TO  LUNAR  EPHEMERIDES 

When  TRACE66  is  operated  in  the  lunar  mode,  it  may  be  necessary  to  cal¬ 
culate  earth-based  observations  of  a  lunar  satellite.  Because  of  interpolation 
errors,  the  order  in  which  such  calculations  are  performed  is  extremely 
important.  Consider,  for  example,  a  lunar  satellite  with  a  selenocentric 
radius  vector  of  2850  km  and  a  geocentric  radius  vector  of  380,  000  km. 

Assume  that  it  is  desired  to  calculate  a  geocentric  range  at  some  time  by 
interpolating  from  an  integrated  moon-centered  trajectory  wich  a  4-min  step 
size  on  the  tabulated  points.  One  possible  approach  would  be  to  transform 
the  moon-centered  coordinates  to  earth-coordinates  using  a  tabulated  lunar 
ephemeris  and  to  interpolate  in  the  earth-centered  system.  Assuming  that 
the  tabulated  epherr.  ris  is  more  accurate  than  the  interpolation.  Fig.  A-l  indi¬ 
cates  that  the  results  should  be  good  to  within  about  10  m.  Alternatively, 
the  interpolation  could  be  done  in  the  moon- centered  system,  before  the 
transformation  to  earth-centered  coordinates  is  performed.  Figure  A-  *  shows 
that  if  the  calculations  are  done  in  this  order,  the  result  is  expected  to  be 
accurate  within  about  0.1  m.  For  the  obvious  reason,  the  latter  method  is 
employed  in  TRACE66. 
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